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ABSTRACT 

Context. Planet traps and snow lines are structures that may promote planetary formation in protoplanetary disks. They are very 
sensitive to the disk density and temperature structure. It is therefore necessary to follow the time evolution of the disk thermal 
structure throughout its viscous spreading. Since the snowlines are thought to generate density and temperature bumps, it is important 
to take into account the disk opacity variations when the various dust elements are sublimated. 

Aims. We track the time evolution of planet traps and snowlines in a viscously evolving protoplanetary disk using an opacity table 
that accounts for the composition of the dust material. 

Methods. We coupled a dynamical and thermodynamical disk model with a temperature-dependent opacity table (that accounts for 
the sublimation of the main dust components) to investigate the formation and evolution of snowlines and planet traps during the first 
million years of disk evolution. 

Results. Starting from a minimum mass solar nebula (MMSN), we find that the disk mid-plane temperature profile shows several 
plateaux (0.1-1 AU wide) at the different sublimation temperatures of the species that make up the dust. For water ice, the correspond¬ 
ing plateau can be larger than 1 AU, which means that this is a snow "region" rather than a snow "line". As a consequence, the surface 
density of solids in the snow region may increase gradually, not abruptly. Several planet traps and desert regions appear naturally as 
a result of abrupt local changes in the temperature and density profiles over the disk lifetime. These structures are mostly located at 
the edges of the temperature plateaux (surrounding the dust sublimation lines) and at the heat-transition barrier where the disk stellar 
heating and viscous heating are of the same magnitude (around 10 AU after 1 Myr). 

Conclusions. Several traps are identified: although a few appear to be transient, most of them slowly migrate along with the heat- 
transition barrier or the dust sublimation lines. These planet traps may temporarily favor the growth of planetary cores. 

Key words. Protoplanetary disks - Planet-disk interactions - Planets and satellites: formation - Planets and satellites: dynamical 
evolution and stability - Accretion, accretion disks - Hydrodynamics 


1. Introduction 

Protoplanetary disk observations provide snapshots of the plan¬ 
etary formation processes that help us understand the physical 
characteristics of such disks. However, understanding planetary 
formation requires modeling the evolution and composition of 
protoplanetary disks. Considering the different physical states of 
the various components of a protoplanetary disk is therefore nec¬ 
essary to constrain favorable scenarios for planetary formation, 
growth, and migration, and solve the apparent inconsistency be¬ 
tween the formation and migration timescales. 

_ Since the discovery of the first exoplanet (Mayor & Queloz 

1 19951) . planetary formation scenarios have f requently been revis¬ 
ited. Even before Charb oimeau et all (I200C ) showed the gaseous 


Goldreich & Tremainel (1 1 9791) . lArtvmowiczl (1 19931) . and Ward 
(1997), who estimated a migration timescale of about 10 5 years 
for a typical Earth-mas s planet in a minimum mass solar nebula 


giant nature of the observed exoplanets, iPollack et alJ (119961) 
described how gaseous planets could form by accretion of 
gas on previously accumulated solid cores of a few times the 
mass of the Earth and estimated a typical planetary forma¬ 
tion timescale around 10 6 ~ 7 years (compatible with the typi¬ 
cal disk lifeti me of a few million years as in ferred from ob¬ 
servati on s by Beckwith & Sargent j 19961) and [Hartmann et al.1 
(Il998l) ). iPanaloizou & Teraueml ( 19991) also noted that hot 
Jupiters were unlikely to have formed in situ, therefore requir¬ 
ing some sort of planetary migration. Type I inward migration 
due to Lindblad resonances with the planet are well known from 


tor a typical Harm-mass planet in a minimum mass solar nebula 
( Weidenschilli ngl fl 9771 : lHavashilll981l) . iKorvcanskv & Pollackl 
(1993) noted that planets should therefore be lost by spiral¬ 
ing into their host star before they could actually grow. Early 
on, IWardl (1199 ll) deta iled the analy t ical e xpression of the coro¬ 
tation torque, before iTanaka et alJ (120021) added a 3D analyti¬ 
cal expression of the vo rtensity gradient h orses hoe drag torque 
(tested numerically by Bate_et al. ([20031) and [D’Angelo et al.l 
(I2003I) ). and lBaruteau & Masset ( 20081) later completed this with 
the entropy gradien t hor seshoe drag torqu e. In the meantime, 
lAlibert et alJ (120051) and llda & Linl (120081) investigated the in¬ 
consistency between the timescales of planet formation and mi¬ 
gration and noted that the mass-distance distribution of the exo¬ 
planets is inconsistent with the rapid inward migration of plan¬ 
etary cores. They noted that the type I inward migration should 
be slowed down by at least an order of magnitude to allow plan¬ 
ets to form and avoid falling into the sta r. Part of the solution, 
investigated by [Hellarv & Nelsonl ( 2012 0. would consist in ac¬ 
celerating the planet formation by considering a proper temper¬ 
ature treatment and a 3D disk model. However, this is not suffi¬ 
cient, and most of the efforts have been concentrated on slowing 
down the inward migration. Various models have been tested to 
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achieve that goal. One of them is the model reported bv lTeraueml 
(200 3]), who studied the effec ts of a regular magnetic field, while 


Nelson & Pana loizou ( 20041) focused on turbulent magnetism. 
Kuchner & Lecar ( 200 2). on the other hand, studied the conse¬ 


quences of an inner cavity and determined that_such a hole could 
halt inward migration. iJang-Condell & Sasselovi (120051) found 
that self-shadowing in the disk can slightly decrease the mi¬ 
gration rate. iMasset et al.1 d2006l) described that a sharp positive 
surface mass density gradient could have similar consequences 
for the creation of planetary tra ps a t the zero-torque radii. Sub¬ 


sequent efforts fro m Paardekooper & Mellema|j20061 2008), 
IPaardekooper & Papaloizoul d2008l). iBaruteau & Masset) (2008), 


IlGev & Cridal d2008l) . iKlev etalJ d2009h . and lAvliffe & Bati 
(1201 Ol) investigated the possibility of slowing down and revers¬ 
ing the inward migration by considering more complete mod¬ 
els involving proper radiative transfer and thermal considera¬ 
tion: this appeared t o be possible for sufficiently low mass plan¬ 
ets (Mp < 40 Mr,). [Menou & Goodmanl (120041) studied the ef- 
fects of opacity transitions on the migration rate and showed 
that if some specific conditions arejnet, the migration could 
he stoppcd in a steady-state a disk. IPaardekooper & Papaloizoul 
(2009a ) estimated that almost any positive surface mass density 
gradient could act as a protoplanet t rap. Zero-torque radii were 
later analyzed by Lyra et a ll (1201 Ol) . IBitsch & Klevl (1201 ll) . and 
lHasegawa & Pudrita (i201 lbl) : they appear to be very important 
in preventing the fall of the planets into their star, but also in cre¬ 
ating a favorable zone for the interactions between planetesimals 
where they are able to combine. iHorn et al.l ( 20121 ) estimated that 
giant planet cores could form at convergent zero-torque radius 
for sub-M 0 planets in 2-3 Myr. 

There appears to be a strong correlation between the migra¬ 
tion rate and the surface mass density and mid-plane tempera¬ 
ture gradients that also strongly depend on the disk composition 
(gas-to-dust ratio, chemical abundances). Although the influence 
of the dust composition is not commonly used in recent numer- 
ical si mula tions of evolving proto planetary disks. Helling et al. 
J20001) and lSemenov et all (12003) studied how opacities are af¬ 
fected by temperature. It therefore appears to be necessary to 
consider how the dust main component phases change in or¬ 
der to estimate the temperature more accurately. From reliable 
evolved disk radial profiles, we skip the planetary core forma¬ 
tion stage and consider how a formed core dynamically inter¬ 
acts with the disk. The resonant torques tha t a planet exerts 
on a disk c an be ca lculated fr om Gold reich & Tremainel (1 1 979t) . 
Ward) (1 1 9881) . and Artvmowiczl (19931, and we considered the 
refinements from IPaardekoo per et al. ( 201 ll) to more accurately 
calculate the various contributions of the corotation torques. In 
their simulations, IBitsch & Kiev! (1201 ll) found a possible equi- 
libriu m radius of a plane t with 2 0 Earth masses around 12.5 
AU. lHasegawa & Pudritzl (1201 lbl) also analyzed the influence 
of the heat transition barrier on the migration torques. Whereas 
most of the previous work has applied the torque formulae 
to simplified semi-analytica l density and temperature profiles 
(lHasegawa & Pudrit3l201 lat IPaardekooper et aUl201 ll) . to sim¬ 
ple densit y prescriptions with a self-consistent 2 D-temperature 
structure (iBitsch & Klevl[20T1] Bitsch et akll2013l) . or to steady- 
state accretion disk models (IBitsch et aTT 2014l)~ we intend to 
apply similar reasoning to more realistic disks obtained from 
numerically simulated viscous evolution rather than analytical 
steady state disks. 

IBaillie & Charnozl (120141) (referred to hereafter as IBC14I) 
have shown that some of the features of viscous a-mode I pro¬ 
toplanetary disks that are observed around forming stars can be 
retrieved numerically using a viscous evolution hydrodynami- 


cal code. These simulations confirmed the importance of jointly 
considering the dynamics, thermodynamics, and geometry of the 
disk. However, as the temperature affects the gas-to-dust ratio of 
the main components of the disk, the opacity of the disk is af¬ 
fected as well. Therefore, it is important to take into account a 
consistent composition of the disk when calculating its tempera¬ 
ture. We here improve the numerical model of BC14 to consider 
these changes in the phases of the disk components. The ther¬ 
mal model includes both viscous heating and irradiation heat¬ 
ing. We follow the evolution of an initial minimum mass solar 
nebula. The geometry (including the delimitations of shadowed 
regions) is calculated self-consistently, with the thermal struc¬ 
ture obtained by semi-analytical radiative transfer calculations. 
The obtained density and temperature radial profiles show dis¬ 
continuities compared to previous results from |BC14| . Density 
bumps mainly result in temperature irregularities such as bumps, 
troughs, and plateaux. The snow or sublimation lines are also 
enlarged. Using our density and temperature profiles, we com¬ 
pute the torque that the disk would exert on a planet embedded 
in the disk. The total torques (resulting from the Lindblad and 
corotation resonances) provide the direction of migration of the 
planetary core within the disk, allowing us to identify conver¬ 
gence and divergence regions, which are also called planetary 
traps and deserts. 

Section |2l details how the hydrodynamical code of 
IBaillie & Charnozj ( 2014 ) is upgraded to consider variations of 
the dust composition with temperature. The density radial pro¬ 
files resulting from the viscous evolution and the calculated ther¬ 
modynamical and geometrical profiles of the disk are shown in 
Sect. [3] These radial profiles, and especially the mid-plane tem¬ 
perature, are analyzed in Sect.[4j where we also discuss the influ¬ 
ence of the composition of the disk and reconsider the definitions 
of snowline and sublimation line. Finally, Sect. [5] calculates the 
resonant torques of potential planetary embryos in the disk and 
follows their migration, defining planetary traps and deserts. 


2. Methods 

2.1. Disk model 

We consider the same model of a v iscous a disk 
(IShakura & Sunvaevl 1 197 3h as was used in IBC14I and use 
the same terminology as in that paper. The turbulent viscosity 
is set to a v isc = 10 2 as is commonly accepted for T Tauri star 
protoplanetary disks without deadzones dFromang & Papaloizoul 
120061) . The time evolution of the surface mass density is given 
by Eq.lTIfrom lLvnden-Bell & Pringleld 19741) . 


dl.(r,t ) 3 d I r d 


dt 


= ~r dr ( ^ dr ^ ^ r ’ ^ 


( 1 ) 


Similarly to BC14, we applied Eq. Q]to a one-dimensional 
grid of masses that are logarithmically distributed in radius be¬ 
tween R * and 1000 AU: each of these masses represents an annu¬ 
lus in the disk. We imposed that the flux at the inner edge cannot 
be directed outward. The inner mass flux gives the mass accre¬ 
tion rate of the disk. The temperature in the mid-plane, T m (r), 
results from the combination of viscous heating, stellar irradia¬ 
tion heating, and radiative cooling in the mid-plane. 

The grazing angle a gr (r ) defines the angle at which the star 
sees the photosphere at a given radial location r. Comparisons 
between an imposed geometry following that prescription and 
a free geometry calculated along with a consistent temperature 
are shown in lBC14l along with a discussion of the necessity of 
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these geometric refinements. The grazing angle, which controls 
the amount of energy that the star provides to the disk photo¬ 
sphere, is related to the photosphere height H p h(r) through Eq. 

El 


a gr (r) = arctan 




- arctan 


H ph (r) - OAR, 


( 2 ) 


A positive grazing angle at a given radius results in an irra¬ 
diated photosphere at that location, while regions that are not 
irradiated are shadowed by inner regions. Using Eq. 18 from 
ICalvet et al.l ( fl99l . we calculated the temperature in the mid¬ 
plane, accounting for viscous heating, stellar irradiation, and ra¬ 
diative cooling. The viscous contribution depends on both the 
surface mass density obtained after temporal evolution and on 
the mid-plane temperature itself (through the viscosity): 


Fy(r) = ±2(r)v(r)^^J = ^l(r)v(r)0 2 (r). (3) 

Therefore, we solved Eq. 18 from ICalvet et al.l (1 199 11 1 nu¬ 
merically as an implicit equation on the mid-plane tempera¬ 
ture. Considering a hydrostatic equilibrium, the vertical density 
distribution follows a G aussian, and we can use Eq. A9 from 
iDullemond et al. ( 200 11 to calculate the ratio x of the photo¬ 
sphere height to the pressure scale height. The disk vertical den¬ 
sity profile is therefore assumed to be the same as an isothermal 
vertical structure at the mid-plane temperature. This approxima¬ 
tion is reasonable below a few pressure scale heights, where most 
of the disk mass is located. The opacities are also functions of 
the temperature, as detailed in Sect. 12.21 We can then estimate 
the corresponding presumed photosphere height H p h at each ra¬ 
dial location and therefore access -jp- Applying Eq. [2] we can 
verify whether the presumed grazing angle has the required pre¬ 
cision or if we should iterate on it. The impossibility of solving 
that problem for any positive value of the grazing angle results in 
a disk column that is not directly irradiated by the star, and there¬ 
fore we removed the irradiation heating term from the mid-plane 
temperature equation. 

Thus, the geometrical structure (photosphere and pressure 
heights) was determined jointly with the temperature by iter¬ 
ating numerically on th e graz ing angle value: the algorithm is 
thoroughly described in lBC141 

2.2. Realistic opacities 

While IbC14I did take into account both sources of heating, they 
only partially considered the thermal variations of the physi¬ 
cal composition of the disk. Indeed, lBC14l considered that the 
dust grains invariably have the same opacity for any tempera¬ 
ture below 1500 Kelvin (considered here to be the sublimation 
temperature of the silicate dust); and another constant opacity 
(a hundred times lower) for temperatures above 1500 K. In the 
present paper, we refine the model by considering a more elab¬ 
orate model of opacities based on optical constants measured 
in laboratory experiments, accounting for the variations of dust 
composition as a function of the local temperature. The com¬ 
putation of t he Rosseland mean opac ities follows the proce dure 
described bv lHelling et al.l<l2000l) and lSemenov et al1(l2003l) . We 
assumed that the dust grains are composed of a mixture of differ¬ 
ent elements: olivine silicate, iron, pyroxene, troilite, refractory 


Elements 

Sublimation 

Temperature 

Relative 

Abundances 

Water ice 

160 K 

59.46 % 

Volatile Organics 

275 K 

5.93% 

Refractory Organics 

425 K 

23.20% 

Troilite (FeS) 

680 K 

1.57% 

Olivine 

1500 K 

7.46% 

Pyroxene 

1500 K 

2.23 % 

Iron 

1500 K 

0.16% 


Table 1 . Sublimation temperatures and relative abundances that affect 
the disk gas opacity. 


and volatile organics, and water ice, with initial abundances (be¬ 
fore volatiles sublimation) given in Table [T] The relative abun¬ 
dances of the opacity model are updated depending on the tem¬ 
perature of the medium and the sublimation temperatures (Table 
[T). We assumed for the dus t grains a modified MRN size distri¬ 
butio n dPollack et alJ 1 19851 IHelling et al.l 12000: Se menov et al.l 
l2003h . that is, grain sizes varying from 0.005 to 5 pm with a 
-3.5 power-law exponent size distribution. The absorption co¬ 
efficients of the composite dust grains were computed using the 
Maxwell-Garnett mixing rule followed by the Mie theory. The 
tabu lated values of the gas opacity computation are taken from 
IHelling et all d2000i) . 


Figure Q] presents the Rosseland (xr) and Planck (k/>) opac¬ 
ity variations with temperature. The Planck mean opacities at 
stellar effective temperature T * = 4000 K in extinction (x P ) 
and absorption ( k * p ) are also shown. These opacities are neces¬ 
sary in p articular for calcula t ing the irradiation heating as de ¬ 
scrib ed in Calvet et al. ( 199il) . lJang-Condell & Sasselovl (l2004l) . 
an d iJang-Condelll (12008" 


It appears that these opacities vary by several orders of mag¬ 
nitude over the concerned temperature range and that these vari¬ 
ations are particularly steep near the sublimation temperature of 
the major disk gas components. These elements and their subli¬ 
mation temperatures are presented in Table [Q sublimation tem¬ 
peratures are given bv lPollack et al.l (1 19941) . corresponding to gas 
densities of about 10~ 1() g.cm~ 3 . 



1 10 100 1000 10000 
Temperature (K) 

Fig. 1 . Mean-opacity variations with local temperature. Black: Rosse¬ 
land mean opacity in extinction. Red: Planck mean opacity in absorp¬ 
tion. Yellow: Planck mean opacity in extinction at stellar irradiation 
temperature. Blue: Planck mean opacity in absorption at stellar irra¬ 
diation temperature. 
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100000 yr 



Disk Midplane Radius (AU) 


(5 points per decade) for an accessible computation time. The 
thermal and geometric structure can then be calculated based on 
the resampled density structure of the disk. 

Assuming a temperature profile following a power law in 
T oc and a temperature resolution of 0.01 K, we can analyt¬ 
ically estimate the highest possible radial resolution for two con¬ 
secutive bins to have temperatures that differ by at least that tem¬ 
perature precision: this optimal resolution is around 100 points 
per decade. To keep a safety margin (in particular to allow for 
the temperature fluctuations due to the opacities), we chose to 
resample our disks with 40 points per decade, which provides a 
good radial precision and allows a good thermal and geometrical 
precision. 


1000000 yr 



Disk Midplane Radius (AU) 


Fig. 2. Compared evolution of the surface-mass density after 100,000 
years (upper panel) and 1 million years of evolution (lower panel) of an 
initial MMSN for various resolutions from 5 to 50 points per decade. 


2.3. Testing the radial resolution 

Our evolution code requires a double check both on the tempo¬ 
ral and on the spatial resolution that the structure of our code 
(density time evolution precedes geometry and temperature cal¬ 
culation and so on for each iteration) allows us to treat sequen¬ 
tially. First, the timestep is adjusted to limit the mass transfers 
to 1% of the available mass in each bin of the simulation, and 
then the radial resolution can be controlled by ensuring that the 
surface mass density profiles are consistent over 1 million years 
for different radial resolutions. Figure [2] shows the radial den¬ 
sity profiles between 1 and 100 AU after 100,000 years and 1 
million years for a range of radial resolutions extending from 5 
points per decade to 50 points per decade. Stronger differences 
around 1 AU are due to the boundary conditions of the test nu¬ 
merical simulations for which the inner edge of the disk was 
considered to be at 1 AU for the purpose of the comparison: as 
the computation time increases exponentially, simulations using 
the highest number of points per decade could only be run over 
a few decades. We note that the surface mass density profiles are 
very similar until 100,000 years, after which the differences are 
more obvious. However, the density profiles remain quite close 
and the only consequence of the difference seems to be a delay 
in the viscous evolution for the simulations with 30 and more 
points per decade. 

Therefore, we estimate that a disk evolution with a good ra¬ 
dial resolution (e.g., 40 points per decade as in the rest of the 
paper) can be approximated from interpolating the density radial 
profile of a disk generated by a less radially resolved simulation 


3. Results 

We foll owed the evolution of an initial minimum m ass solar 
nebula (IWeiden schilling 119771) with the scaling from Havashi 
(119811) : 


£(r) = 17,000 (y^j) 3/2 kg ■ nr 2 . (4) 

The central star was a classical T Tauri type young star 
with constant M* = 1M Q , R , = 3 R Q , 71, = 4000 K and 
X* = 4n 7? 2 ct b 7ij throughout the simulation. 

The choice of the minimum mass solar nebula (MMSN) is 
motivated by several arguments. First, 1BC14I investigated a di¬ 
versity of initial disk conditions (total disk mass and radial dis¬ 
tribution of the initial angular momentum) and found that these 
numerical simulations converged to similar steady states (char¬ 
acterized by a uniform mass flux and an asymptotic surface 
mass density distributio n in r -1 ), although in sli ghtly different 
timescales. In addition, IVorobvov & Basul ( 2007 ) showed that 
the MMSN density profile (following a power-law in r 15 ) was 
consistent with an intermediate stage of an evolving protoplane¬ 
tary disk under self-regulated gravitational accretion. Therefore, 
the MMSN profile makes an initial profile as reasonable as any 
other snapshot that could have been taken in the disk evolution. 
Finally, this fiducial case makes sense because we can better 
compare our results with previous studies. 


3. 1. Time evolution 


Figure [3 presents the evolution of such a disk over 10 mil¬ 
lion years. Although the protoplanetary disk gas is believed 
to photo-evaporate in a few million years (Font et all [2004; 
lAlexander & Armitagel 120071120091: lOwen et al.11201 0l) .~ we let 
our numerical simulations extend over longer times to reach a 
stationary steady-state, which is characterized by a uniform mass 
flux (within one order of magnitude maximum), as seen in Fig. 

SI 

As IBC14I detailed, the disk starts to spread viscously out¬ 
ward before beginning to accrete onto the central star after a 
few thousand years. Figure [3 shows that the surface mass den¬ 
sity radial profile can be modeled as a power law by segments. 
While two segments seem to be sufficient until 10,000 years, 
longer evolution times require at least one more segment. In the 
early times, the single connexion point between the two power- 
laws connects an inner region that has already evolved toward 
a steady-state power-law radial profile for this region, and an 
outer region that is much closer to the initial state. We therefore 
call that radial location the relaxation radius. After 10,000 years. 
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0.1 1 10 100 1000 
Disk Midplane Radius (AU) 


Fig. 3. Surface mass density radial profile evolution for an initial min¬ 
imum mass solar nebula in the case of a self-consistently calculated 
geometry with a full continuous model of opacities. Simulations are re¬ 
sampled to 40 points per decade. 


the (inner) relaxation radius is located between 7 and 10 AU. 
After 1 million years, a secondary "knee" appears around 250 
AU. In the intermediate region (10-250 AU), the viscous evo¬ 
lution appears weaker than in the simulations of 1BC14I with a 
simpler opacity model. At 1 Myr, the surface mass density pro¬ 
file can be modeled as a power law between 10 and 250 AU 
with E(r) oc r l2 , and at 10 Myr, we can approximate the den¬ 
sity profile by X(r) oc r 11 . In this region, and evolved disks, 
the power-law indices are comparable with those from 1BC14I 
and with t hose observed by Isella et al. ( 20091) in the Tau rus re¬ 
gion and lAndrews et alJ (l2009h and Andrews et al.1 d2010l) in the 
Ophiuchus region. 

Figure [4] shows the time evolution of the mass flux as a func¬ 
tion of the radial distance. We note thatjhese radial profiles 
present bumps that were not visible in IBC141 Although the di¬ 
rections and amplitudes remain globally identical, some irregu¬ 
larities appear; for example, in the 7-10 AU region where the 1 
Myr flux profile is briefly inverted and the 10 Myr profile is no 
longer uniform over that range (the flux varies by almost a factor 
10 ). 



0.1 1 10 100 1000 
Disk Midplane Radius (AU) 

Fig. 4. Mass flux radial profile evolution for an initial minimum mass 
solar nebula in the case of a self-consistently calculated geometry with 
a full continuous model of opacities. 


3.2. Thermal evolution 

The thermal evolution of the disk is presented in Fig. [5] The 
mid-plane temperature is calculated in a similar way to lBC14t 
However, the calculation of the mid-plane temperature takes into 
account a composition of the dust that is consistent with the tem¬ 
perature (i.e., where sublimated species have been removed from 
the opacity calculation according to TableQ} while iterating over 
the possible grazing angles and mid-plane temperatures: the tem¬ 
perature obtained after the algorithm converged therefore vali¬ 
dates the consistent composition. As described in Sect. [2] the 
geometric and thermal structure are recalculated over 40 points 
per decade. 

As in IBC141 the mid-plane temperature presents irregular 
features that did not appear in the radial profiles obtained with 
the simpler opacity model of BC14L In particular, we note the 
temperature plateaux at the change of phase temperatures of the 
dust components (see Table Q}. As the disk evolves, the sub¬ 
limation plateaux drift inward: the silicate sublimation plateau 
is inside 0.2 AU after 1 Myr while it was originally between 
0.3 and 1.3 AU. In addition to these plateaux, we note the 
troughs around 10 AU and 250 AU, which drift inward as the 
disk evolves. These features are analyzed more thoroughly in 
Sect. 14.21 However, it appears that the temperature in the re¬ 
gion between these features can be modeled as a power law, with 
T(r) oc r" 0 47 beyond 1 Myr, which rec alls the usual app r oxima - 
tion of IChiang & Goldreichl d 19971) or iDullemond et akl d2001l) . 
according to which T(r) oc r - ^ 2 . 



Disk Midplane Radius (AU) 

Fig. 5. Evolution of the mid-plane temperature as a function of the ra¬ 
dial distance for an initial minimum mass solar nebula in the case of a 
self-consistently calculated geometry with a full continuous model of 
opacities. 


3.3. Geometry evolution 

As the geometry of the disk is calculated jointly and self- 
consistently with the thermodynamical structure, irregular fea¬ 
tures can also be expected in the height radial profiles, as shown 
in Fig. [6j the pressure scale height shows radial gradient dis¬ 
continuities similar to those observed in Fig. 0 Similarly, it is 
possible to approximate the pressure scale height by an almost 
constant power law between 20 AU and 250 AU: h piessme = r 126 
after 1 Myr and /r pr essure = r 1 ' 28 at 10 Myr. These index values are 
very c lose to the 9/7 index suggested bv IChiang & Goldreichl 
( 1997 ) in the case of a passive disk for which the photosphere 
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height (height at which the line-of-sight optical depth equals 1) 
would be directly proportional to the pressure scale height. It ap¬ 
pears then that the approximation of h pressur:e oc r’ n can only be 
valid in the region between 10 and 250 AU, where the disk pho¬ 
tosphere is irradiated and the opacity varies smoothly with the 
temperature (i.e., where the temperature is lower than the subli¬ 
mation temperature of the main components, as listed in Table 

m. 


iDullemond et al. ( 2001 1. the values are in the range from 1.5 to 
6. However, unlike the usual approximation (x = 4) suggested 
by lChiang & Goldreichld 1997b . this ratio is neither constant nor 
uniform. This x profile, like the grazing angle a profile (Fig. H 
show similar gradient discontinuities as in the surface mass den¬ 
sity or temperature radial profiles. The non-irradiated zones also 
appear quite obvious in these figures because they correspond to 
locations where ^ and the grazing angle are not defined. 



0.1 1 10 100 1000 
Disk Midplane Radius (AU) 



Disk Midplane Radius (AU) 


Fig. 6. Evolution of the pressure scale height as a function of the ra¬ 
dial distance for an initial minimum mass solar nebula in the case of a 
self-consistently calculated geometry with a full continuous model of 
opacities. 


The radial profile of the photosphere height (Fig. [7} reveals 
a zone inside of 10 AU where the disk photosphere is irregu¬ 
larly irradiated, while it is entirely irradiated in the outer regions. 
Between 20 and 250 AU, the pressure scale height can be ap¬ 
proximated by a power law: // p i loto oc r u 1 after 100,000 yr and 
^pressure x r 114 at 1 Myr and 10 Myr, These index values are 
reminiscent of those expected by [Kenyon & Hartmann! (1 987 ) 
around 9/8 and obtained by the numerical simulations of BC14 . 
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Fig. 8. Evolution of the pressure scale height to photosphere height ratio 
as a function of the radial distance for an initial minimum mass solar 
nebula in the case of a self-consistently calculated geometry with a full 
continuous model of opacities. Missing points are due to the regions of 
the disk that are not directly in the stellar line of sight at a given location 
and evolution time. 



Disk Midplane Radius (AU) 

Fig. 9. Evolution of the grazing angle as a function of the radial dis¬ 
tance for an initial minimum mass solar nebula in the case of a self- 
consistently calculated geometry with a full continuous model of opaci¬ 
ties. Missing points are due to the regions of the disk that are not directly 
in the stellar line of sight at a given location and evolution time. 


Fig. 7. Evolution of the photosphere height as a function of the radial 
distance for an initial minimum mass solar nebula in the case of a self- 
consistently calculated geometry with a full continuous model of opac¬ 
ities. 

The evolution of the photosphere height over pressure 
scale height ratio x is displayed in Fig. [8] As predicted by 


4. Discussion 

4. 1. Influence of the disk composition on density and 
temperature local gradients 

By comparing these results with the simulations of IBC141 we 
note how important it is to take the physical composition of the 
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dust into account by using a complete model of opacities that 
varies with temperature. The main difference in the surface mass 
density radial profiles resides in the secondary knee th at appe ars 
at 1 Myr around 250 AU. Earlier conclusions from IBC14I re¬ 
main valid: the surface-mass density profile evolves toward a 
shallower profile, first in the inner regions and then increasingly 
farther away from the star. However, when the relaxation radius 
reaches 250 AU, the profile is then fragmented into three power- 
law segments. This creates a discontinuity in the density gradient 
that may strongly influence the exchange of angular momentum 
with possible planetary cores within the disk (see Sect. |5j. In 
addition to that, we note that the mass flux may reverse in the 
middle of the disk until as late as 1 Myr. These density bumps 
obviously have consequences for the thermodynamical and geo¬ 
metrical structures. The mid-plane temperature and grazing an¬ 
gle radial profiles do in fact present new features, particularly 
at the radii of these discontinuities, when the physical compo¬ 
sition of the disk is properly treated. Temperature plateaux are 
other important consequences as they may affect, for instance, 
the snowline position (see Sect. ro . 

4.2. Temperature radial profile analysis 

Figure [10] details the temperature structure in the mid-plane of 
the protoplanetary disk after 1 Myr of evolution. As well as the 
mid-plane temperature. Fig. [TO] also displays the grazing angle 
and optical depth radial profiles, and the distribution of the heat¬ 
ing sources. 


1 Myr 



0.1 1 10 100 1000 
Disk Midplane Radius (AU) 


Fig. 10. Mid-plane temperature radial profile (black) after 1 mil¬ 
lion years of evolution of a minimum mass solar nebula with a self- 
consistently calculated geometry and a full continuous model of opaci¬ 
ties. Disk shadowed regions are displayed in gray. The ratio of the vis¬ 
cous heating contribution over the total heating (viscous heating rate) 
is presented in red, the grazing angle radial profile in yellow, and the 
optical depth radial profile in blue. 

We first note is that the black curve shows temperature 
plateaux that coincide with the sublimation temperatures listed 
in Table Q] These plateaux are not flat, but present variations of 
a few Kelvin over radial regions that may be up to 1 AU wide. 
As Fig. [T] showed, opacities may vary quite strongly with tem¬ 
perature around these changes of phases. Thus, the temperature 
drop induced by a local surface mass density variation can be 
compensated for by an opacity variation over a few Kelvin to 
provide the equivalent heat and maintain a quasi-constant tem¬ 
perature. In these regions, the irradiation heating effectively has 
a more important contribution than outside, where it tends to be¬ 
come negligible compared to the viscous heating. This is con¬ 
firmed by the red curve, which represents the ratio of the viscous 
contribution over the total heating received by the disk (viscous 


heating rate): the troughs at 0.1, 0.4, 0.7, 1.2, and 2 AU coincide 
with the tempera ture plateaux. This transition i s called heat tran¬ 
sition barrier by [Hasegawa & Pudritz] (1201 lbl) . In addition, we 
note that the silicate sublimation plateau migrates inward as the 
disk evolves, and after 1 Myr, it is located at the disk inner edge. 
We interpret the inner disk physical structure as follows: 

- The surface mass density is generally a decreasing function 
of the radial distance. 

- When a phase-transition temperature is reached, the opacity 
suddenly increases because of new condensed species, even 
though the temperature remains stable over a few fractions 
of an AU. 

- The viscosity v = o\,i sc c s /z pressure therefore follows the tem¬ 
perature stability. As the density and the angular velocity de¬ 
crease outward, the quantity of viscous heating decreases as 
well. 

- To maintain the temperature, the irradiation heating tries to 
compensate for the viscous heating loss, looking for an opti¬ 
mal grazing angle that can maximize the irradiation heating 
from the star. 

- At some radii, the disk photosphere geometry is such that the 
irradiation heating can no longer compensate for the loss in 
viscous heating. The grazing angle then drops until the disk 
photosphere is not irradiated anymore. The viscous heating 
remains the only source of heating, resulting in a shadowed 
region. 

- The temperature decreases again because of the lack of irra¬ 
diation heating. However, the opacity variation is smoother 
now that we are no longer at the phase transition temper¬ 
ature. The viscosity now decreases, reinforcing the decrease 
in viscous heating. In the meantime, the opacity is now much 
higher than it was before the change of phase. 

In addition to these plateaux, the mid-plane temperature ra¬ 
dial profile shows two important troughs with an amplitude drop 
of about 10 Kelvin. The first one, located around 10 AU, co¬ 
incides with the limit where the dominating source of heating 
changes: viscous heating is clearly stronger below 10 AU, while 
stellar irradiation heating dominates outside that limit (red curve 
from Fig. ITOl). The second trough, around 250 AU, generates a 
drop in temperature from 20 to 8 Kelvin. At these temperatures, 
the opacity varies by two orders of magnitude over a temperature 
range of 20 Kelvin (Fig. [T]): the grazing angle drops and the ge¬ 
ometry tends to shadowing. The temperature bump disappears in 
simpler opacity models (IBC14I) . In addition, it remains present 
(although less pronounced) in earlier evolution ages. 

Furthermore, we note that in the irradiation-dominated re¬ 
gion, the disk is permanently irradiated, which is not the case in 
the inner regions (see gray bands in Fig. ITOt . These shadowed 
regions coincide with a drop in the grazing angle (see yellow 
curve in Fig. [13. Moreover, all the temperature plateaux (all are 
located below 10 AU) correspond to irradiated zones, while the 
outer edges of the plateaux (where the temperature starts drop¬ 
ping again) trigger shadowed regions, as detailed above. The 
regio n between 20 and 250 A U matches the two-layer model 
of IChiang & Goldreichl d 1997b well because the temperature is 
clearly below the dust sublimation temperatures, and therefore 
in a temperature domain where the opacity varies so smoothly 
that it can approximated by a constant, allowing the retrieval of 
the simulation results from lBC14l 

Finally, the stellar line-of-sight optical depth radial profile 
(blue curve in Fig. ITOl) shows that the disk becomes optically 
thin beyond 20 AU. However, the optical thickness drops again 
at the heat transition barrier (around 10 AU at 1 Myr). 
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4.3. Snow region 

While the snowline is defined as the radius at which the tem¬ 
perature is equal to the water ice condensation temperature (160 
K), the fact that we now observe a plateau around that tempera¬ 
ture shows that we should rather talk about a snow region than a 
snowline. Figure|TT]presents the evolution of the inner and outer 
edges of the plateau around the water ice condensation temper¬ 
ature. Given the precision on the temperature, and in the opac¬ 
ity model, we chose to define the plateau temperature range as 
160 K ± 2 K. The snow region can be as wide as 1 AU at early 
ages and migrates inward until it stabilizes below 2 AU in a few 
Myr. 


6 


13 s 

CD 

c 4 

o ^ 

N 

<D 

u 3 

(D 

1 2 

1 

10 1 10 2 10 3 10 4 10 5 10 6 10 7 
Evolution Time (yr) 

Fig. 11. Time evolution of the snow region (mid-plane radial location 
for which the temperature coincides with the water-ice condensation 
temperature ± 2 K). 

Similarly, we define the silicate sublimation region in place 
of the silicate sublimation line at 1500 K + 20 K, and we present 
the evolution of its edges in Fig. [12] The sublimation region ini¬ 
tially spreads from 0.3 to 1.3 AU. It also migrates inward and 
seems to reach the inner edge of the disk in a few Myr. 




Fig. 12. Time evolution of the silicates sublimation zone (mid-plane ra¬ 
dial location for which the temperature coincides with the silicate sub¬ 
limation temperature ± 20 K). 

It appears that the zone of the mid-plane in which the sub¬ 
limation of an element can occur is AU-wide instead of the ex¬ 
pected sharp frontier usually called snow or sublimation line. 


At this temperature, the corresponding element is sublimated, 
which affects the medium opacity and therefore the heating re¬ 
ceived by the disk at this location. This heating variation is 
compensated for by a grazing angle variation that maintains a 
smoother continuity of the temperature and gas-to-dust ratio. 
The origin of the plateau is therefore related to the partial and 
gradual sublimation of the element that is sublimated at this tem¬ 
perature: the gas-to-dust ratio of that element is close to 0 at the 
outer edge of the plateau while it increases inward up to 1 at 
the inner edge. The other elements gas-to-dust ratios are mostly 
unaffected at this temperature. Partial sublimation could occur 
layer by layer on dust particles or on a vertical scale (the column 
height of the sublimated material increasing inward). 

The changes of phase can now occur in much wider regions 
than the previous snowlines. This may favor smooth variations 
in physical compositions across the disk as the planetary embryo 
migrates inward across the snow or sublimation region. This may 
influence the chemical models, for example, those that try to ex- 
plai n the abundan c e of ca rbon in the so-called carb o n-rich plan¬ 
ets: iFortnev et al] ([2010) and Madhusudh an et alJ (1201 ll) sug¬ 
gested that, at equilibr ium, all the available O should go into or¬ 
ganics, whereas iVenot et all ( 2012 ) described a non-equilibrium 
chemistry to model the C/O ratio. Therefore, the shallow varia¬ 
tions of the temperature profile may favor equilibrium chemistry. 


5. Consequences on planet traps and deserts 

The results reported above clearly show that rapid variations of 
temperature and density occur in the disk in the transition region. 
It is thought that such a transition may potentially create planets 
traps and deserts. This has been explored previously in static- 
disk models, but never in dynamically evolving disks. Thus in 
the following we compute the torque that the disk would exert 
on a putative planet. In particular we explore how planet traps 
appear, move, and disappear in the disk. 

Assuming that it has already formed, a planetary 
embryo exchanges a n gular momen tu m w it h the disk 
J Goldreich_&Tremaine _ J_979[_ JWarcj 119881 : lArtvmowiczl 


119931 : Jang-Condell & Sasselovl l2Q05t> . These exchanges are 


due to the resonances excited by the planetary embryo in the 
disk (Lindblad resonances caused by the action of the induced 
spiral arms, and corotation resonances). Thus, the planet exerts 
a torque on the disk and therefore the disk exerts an opposite 
torque on the planet. We calculate these torques in the case 
of the evolved disk described in the previous section. We then 
study their effects on potential planetary embryos. We assume 
here that the disk structure is not modified by the planet. 


5.1. Lindblad torques 

A given perturber, such as the planet in our disk, excites 
Lindblad resonances of multiple orders iGoldreich & Tremaine! 
Il979[ il98 0). Using a two-dimensional approximation, con¬ 
sidering laminar disks, a planet on a circular orbit, ignor¬ 
ing the disk self-gravity and assuming thermal equilibrium, 
[Paardekooper & Panaloizoul ( 2008 ’) were able to derive the fol¬ 
lowing formula for the total Lindblad torque exerted by the disk 
over the planet: 


1 'l.indblad — 


r 0 (rp) 

r 



- 1.7 


dlnT 
<3 In r 


(rp) + 0.1 


d In 2 
<91nr 



with y = 1.4, the adiabatic index, 
r 0 (rp) = (f) 2 20>)r 4 (Q0>)) 2 , 


(5) 
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^ _ Vess(rp) 

~ r P 

and Q.(rp) the Keplerian angular velocity at the planet position 
in the disk. 

Although we just showed in the previous sections that the 
disk can present radial variation in density and temperature gra¬ 
dients, the tor que expressi o ns are still evaluated at the pla net ra¬ 
dial location. Wardl (119971) . iHaseeawa & Pudrit3 (1201 1 alibi) , and 
iMassetl (1201 1 ) (Appendix B) developed fully analytic torque ex¬ 
pressions that could account for the local gradients at each res¬ 
onance location rather than just considering the gradients at the 
planet radius. However, these formulas could not retrieve the am¬ 
plitude of the expression derived in Paardekooper & Papaloizoul 
( 120081 ) . which was benchmarked against numerical simulations. 
In addition, more subtle expressions involving the second deriva¬ 
tive of the temperature derived by iMassel ( 201 111 (Eq. 79) were 
only tested in an isothermal disk. To be consistent wi t h the coro- 
tation torque expressions from iPaardekooper et akl (1201 ll) . we 
therefore estimated the total Lindblad torque using the expres¬ 
sion of iPaardekooper & Papaloizoul ( j2008i) (our Eq. |5). 

We note that the Lindblad torque presents a stronger depen¬ 
dence on the temperature gradient than on the density gradient, 
which reflects the effect of the pressure buff er described i n Ward 
(119971). Following a similar pr ocess as in Bitsc h & Klevl (2011) 
and Bitsc h et alJ d2013l 120141) . we can use here the results of 
Sect. [3 which provide the density and temperature of an evolved 
disk at a given date, therefore gaining in consistency as the tem¬ 
perature is not set to a power law and is calculated jointly with 
the geometry of the disk from the density resulting of the viscous 
evolution. 


5.2. Corotation torques 


Corotation resonances are known to exert compli¬ 
cated torques_ tha_t_ include linear and nonlinear parts. 
iPaardekooper & Papaloizoul ( 2009bl) showed that the coro¬ 
tation torques are generally nonlinear in the usual range of 
viscosity (cr,j SC < 0.1). The nonlinear contribution, due to the 
horseshoe drag dWardlll99H) caused by the interaction between 
the planet and the fluid element moving in its vicinity, is also 
known for having two possible origins: barotropic, intially 
formalized by T anaka et akl (120021) . and entropic, detailed by 
iBaruteau & Masset (2008). 

Concerning the horseshoe drag, IPaardekooper et akl (1201 ll) 
described the density perturbation generated by the corota¬ 
tion resonances and provided expressions for both the entropy 
and vortensity (or barotropic) contributions. Assuming a grav¬ 
itational softening b = 0.4/i prcss , iBitsch & Klevl (1201 ll) and 
Bitsc h et aD (120141) summarized these expressions to obtain the 
following contributing torques: 


It appears that this unsaturated corotation torque strongly 
depends on the temperature and surface mass density gradi¬ 
ents. It also scales with /Vfj,_as does the Lindblad torque. How¬ 
ever, iPaardekooper & Papaloizou (2009a) showed that given 
the viscous, diffusive, and libration timescales, the linear ef¬ 
fects of the corotation torques can be saturated for some vis¬ 
cosities and some planet masses. For our disk that evolved 
for 1 Myr, the viscosity range compared to Fig. 14 from 
iPaardekooper & Papaloizoul d2009al) suggests that saturation 
cannot be neglected for planetary masses higher than 6 M e . 
IPaardekooper et akl (1201 ll) defined weight functions for the par¬ 
tial saturation of the corotation torque. These functions vary with 
the half-width of the horseshoe, which depends on the mass of 
the planet. Appendix A of iBitsch & Kiev! (1201 ll) summarized 
this method and added correcting factors. We used a similar 
torque calculation, which is necessary to take into account the 
variations with the planet mass. 

5.3. Planetary core migration 
5.3.1. Migration direction 

In the present exercise, we set the planet mass to a typical giant 
planet core mass Mp = 10M© and varied rp, the initial planet 
location. The total torque exerted by the disk over the planet is 
given by 


r\ot — ILindblad T Tcorotation- (9) 

A positive torque exerted by the disk on the planet means that 
the disk gives angular momentum to the planet and therefore that 
the planet migrates outward. In contrast, a negative total torque 
results in the disk gaining angular moment from the planet, and 
the planet migrating inward. Thus, at the locations where the 
total torque sign changes, we can define two physical radii: 

- If the total torque is negative inside and positive outside that 
line, the potential planets are driven away from that loca¬ 
tion, therefore defining a depleted region in planetary cores, 
corres ponding to the planet deserts of Has egawa & Pudrita 
( 120121 ) . 

- If the total torque is positive inside and negative outside, 
the potential planets converge toward this location, which 
we thus call a planetary trap. Such a zero-torque radius 
was called bv lLv raetal.1 ( 20101) the equilibrium radius. 
[Hasegawa & Pudritzi (1201 la ) studied the vertical effects of 
the planet mass, but only considered fixed surface mass den¬ 
sity radial profiles, while lHasegawa & Pudritz ( 1201 lb ) esti¬ 
mated planet trap locations due to Lindblad torques alone. 


r hs,entro — 

_ r oO>) 7 9 

r 

r hs,baro — 

- r ^i.i| 


dlnT <91nl 

.( rp ) + ( r -l)__( rp) | (6 ) 


dim 


dim 


(7) 


iPaardekooper et al . ( 201 1) showed that in the absence of sat¬ 
uration of the linear contributions, the total corotation torque can 
be defined as 


Lcomtation — I hs.haro T I hs.entro- 


( 8 ) 


5.3.2. Density and temperature features at the origins of the 
traps and deserts 

Figure [13] presents the radial profiles of the Lindblad, corota¬ 
tion and total torques after 1 million years of evolution, when 
the planetary cores can already exist but the gas disk is not yet 
completely photo-evaporated. We chose to focus on the plane¬ 
tary formation region and therefore to limit the radial extent of 
our migration investigation up to 20 AU. The torques are nor¬ 
malized by To. A left-pointing arrow shows a negative torque 
and therefore an inward migration, while a right-pointing arrow 
shows a positive torque or outward migration. 

The Lindblad torque is mostly negative except around 10 
AU, where the temperature gradient is reversed. The corotation 
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cores at radii compatible with known gas giant planet radii. This 
is compatible with our present study and the mass accretion rates 
obtained by |BC14l for evolution ages younger than 1 Myr. 
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Fig. 13. Radial radial profile of the migration torques (black arrows 
showing the direction of migration: outward for positive torques and 
inward for negative torques) after 1 Myr of evolution. Lindblad torques 
(blue), corotation torques (red), and shadowed regions (gray) are also 
represented. 


Fig. 14. Radial profile of the total torque (black arrows showing the di¬ 
rection of migration) after 1 Myr of evolution. Temperature radial pro¬ 
file (blue curve), surface mass density radial profile (red curve), and 
shadowed regions (gray regions) are also represented. 


torque is mainly positive, but it can become negative in the tem¬ 
perature plateaux regions where the temperature gradient is shal¬ 
lower (Fig. 1141). The orders of magnitude of the two contributions 
are comparable, and the resulting total torque alternates between 
positive (outward) regions and negative (inward) regions. After 1 
million years of evolution, outward migration could occur below 
0.33 AU, between 0.41 and 0.60 AU, between 0.85 and 1.7 AU, 
and between 2.2 and 8.0 AU. As a consequence, planets could 
accumulate at the traps located around 0.33, 0.60, 1.7, and 8.0 
AU and the regions around 0.41, 0.85, and 2.2 AU should be 
strongly depleted in planets. The snowline (dotted line at 1.8 AU 
in Fig. [T3l > seems to be closely related to the nearby trap at 1.7 
AU: the inner edge of the 160 K-plateau coincides with a trap, 
while there is a planetary desert at the outer edge of the plateau, 
as shown in Fig.[l4j which presents the total torque together with 
the surface mass density and temperature radial profiles. As a 
comparison. iBitsch & Klevl (1201 li) found a possible equilibrium 
radius around 12.5 AU for a 20 Earth-masses planet, which could 
be consistent with our 8 AU-trap. 

As the total torque strongly depends on the surface mass den¬ 
sity and temperature gradients (see Eqs. m , we can here relate 
the torque radial profile with the bumps in surface mass density 
and troughs in temperature. We note here that the density bump 
around 6 AU seems to result in a stronger posi tive torque gradi¬ 
ent, in agreement with the estimate of iMasset et ali ( 200 6) that 
sharp surface mass density gradients are required to slow-down 
type I migration. On the other hand, the temperature irregular¬ 
ity located between 9 and 10 AU seems to be at the origin of 
the torque inversion in this region. The density profile there does 
not seem to affect the total torque trend. This sharp tempera¬ 
ture gradient appears to be correlated with an inner shadowed 
region, as expected by iKretke & Linl ( 2012 b. who described a 
model in which the absence of stellar irradiation could lead to 
an outward migration, and therefore estimated that sustaining an 
outward migration in an irradiated and active disk would be dif¬ 
ficult without any shadowing effects. This correlation confirms 
the importance of properly taking into account the geometry and 
photosphere irradiation. [Kretke & Linl (12012|) found traps inside 
1 AU for the most massive embryos and estimated that a mass 
accretion rate M > 10 s M 0 .yr” 1 was required to trap planetary 


5.3.3. Earlier time evolution 

To follow the evolution of the trap and desert locations. Fig. ED 
presents the total torque profiles as functions of the radial dis¬ 
tance for different ages. After 1000 years, we can only find two 
traps around 0.22 and 11 AU and one desert located at 7 AU. The 
trap is correlated with the shadow region just outside the temper¬ 
ature plateau at the water ice sublimation temperature (160 K). 

After 10,000 years, we now have two traps (4 and 13 AU) 
and three deserts (3.5 and 6 AU), while the water ice sublimation 
plateau appears to be located around the 5 AU. The outer trap 
may be due to the temperature gradient inversion there. 

After 100,000 years of evolution, we now have four traps 
(0.45, 0.75, 2, and 10 AU) and four deserts (0.23, 0.6, 1, and 3 
AU). The traps appear to be located in the outer parts of the shad¬ 
owed regions, while three of the deserts are at the inner edges 
of the shadowed zones. In addition, the water ice line coincides 
with the trap at 2 AU. 

Traps and deserts seem to relate better to the variation in the 
temperature gradient than to the surface mass density gradient 
trends, as expected from the coefficients in Eq.0 


5.3.4. Trap and desert migration 


To better follow the evolution of the traps and deserts. Fig. [16] 
presents the variations in number and position of these traps and 
deserts. 

Figure IT6l shows the evolution of the traps and deserts. We 
first note that in the early evolution of the disk, it seems to be 
possible to create traps and deserts below 0.3 AU. However, as 
these traps disappear in a few thousand years (much before a 
steady-state is reached), they probably only exist as long as the 
disk "remembers" the initial condition of the simulation. Sec¬ 
ond, it appears that there is a permanent trap correlated with the 
heat transition barrier (dashed line): this trap is systematically 
slightly (less than 1 AU) inward of that viscous-irradiation tran¬ 
sition (consistent with the planetary trap p osition of a 20 Earth- 
masses pla net estimated by Bitsch & Klevl (1201 lh ). recalling the 
estimate of lHasegawa & Pudritzl (1201 lbb that the heat transition 
barrier was at the origin of trapping regions. Between 10,000 
and 100,000 years, the heat transition barrier moves slightly out- 
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Fig. 15. Radial profiles of the total torque (black arrows showing the 
direction of migration) after 1000 years (upper panel), 10,000 years 
(middle panel), and 100,000 years (lower panel) of evolution. Tempera¬ 
ture radial profiles (blue curve), surface mass density radial profiles (red 
curve), and shadowed regions (gray regions) are also represented. 


ation radius that crosses the planetary formation region at this 
age (see Sect. |3}. That trapping location is the outer bound¬ 
ary of an outward-migrating region that initially spreads from 
around 6.5 AU (about 1 AU outer to the water ice sublima¬ 
tion line) up to ~ 13 AU. After 10,000 years of evolution, that 
outward-migrating region is divided into two (sometimes three) 
sequences of (outward-migrating + inward-migrating) regions. 
A depleted zone in planets can be found ~ 1 AU outside of the 
water ice sublimation line at any time. 

After 10,000 years, we also note traps and deserts following 
the main dust component sublimation lines: 

- the water ice sublimation line coincides with a systematic 
trap, 

- the volatile organics sublimation line is closely accompanied 
by an inner trap and an outer desert, 

- the refractory organics and troilite sublimation lines are sim¬ 
ilarly bordered, if not as closely as the volatile organics line, 

- the silicate sublimation line is closely followed by an outer 
desert region. 

From Fig. [14] we can infer that the space between two sub¬ 
limation regions is a zone of outward migration: therefore the 
chances of a planet to become trapped or significantly slowed 
down in the inner disk are not negligible. Once trapped, the form¬ 
ing planet is very likely to follow the trap migration. In addi¬ 
tion, most of the planet traps (those that are correlated with the 
heat-transition barrier or the ice lines) seem quite sustainable: 
these traps and deserts slowly migrate inward with the sublima¬ 
tion lines as the disk ages and cools down. Only a few isolated 
transient traps and deserts appear between the water-ice line and 
the heat-transition barrier. However, it is very likely that a planet 
in such a disappearing trap would be trapped again in the next 
outer trap. 

iKlev et all ( 2009 ) suggested that a possible scenario to ad¬ 
dress the inconsistency between the formation and migration 
timescales could be the retention of the icy cores at the snow¬ 
line. We extend this: it is possible to trap planets at any main 
dust component sublimation line. Dust sublimation is a physi¬ 
cal mechanism sufficiently efficient to trap planets along the ice 
lines. In addition, as suggested by the correlation of the trap 
or desert locations with respect to the shadowed regions (Fig. 
M- the heating sources play an important role in trapping plan¬ 
ets. This is reinforced by the presence of traps along the heat- 
transition barrier. 



Time (yr) 


Fig. 16. Time evolution of the migration traps (blue "+") and deserts 
(red "x"). The snowline position (dotted line) and the heat transition 
radius (dashed line) are also represented. 


ward, which may coincide with the outward drift of the relax- 


5.3.5. Perspectives 

Although multiple planet-systems seem to be the most common 
configuration, we here only studied the interaction of a single 
planet with the disk. ICossou et al. ( 2013 ) showed that multiple 
planetary cores tend to be trapped in chains of mean motion res¬ 
onances. Therefore, a proper treatment of multiple planet inter¬ 
action with the disk requires considering the dynamics of these 
planets in addition to the disk dynamics. Coupling our hydrody- 
namical approach with an N-body code requires a more thorough 
study and will be the object of a future paper. In addition, as we 
have detailed, our current model does not take into account the 
feedback of the planet on the disk, which is currently not affected 
by the planet. This might also be necessary for future studies. 

We leave the consideration of dead zones for a future paper 
because they probably generate density and temperature irregu¬ 
larities that have been shown to be favorable for creating plane¬ 
tary traps and deserts. Although the torque expressions detailed 
above show a limited dependency on the planet mass (oc /V/p,)), 
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this is kept for a subsequent review, along with the dependency 
on the disk mass and the mass accretion rate. 

6. Conclusions 

We have studied the viscous evolution of a protoplanetary disk, 
using an a prescription for its viscosity. Starting with an ini¬ 
tial minimum mass solar nebula, we self-consistently calculated 
the geometry, thermodynamics, and density distribution evolu¬ 
tion of such a disk, taking into account the various sources of 
heating (viscous and irradiation), the shadowing effects, and the 
physical composition. The physical phases of the various dust 
components have a dramatic influence on the opacity of the disk 
and therefore reflect on the mid-plane temperature. The evolved 
density radial profiles present a few irregularities compared to 
the previous results of iBC14l These density bumps mainly re¬ 
sult in temperature features (bumps, troughs, plateaux) and also 
widen the snowline and sublimation line regions. Our main re¬ 
sult is that snowlines may better be called snow regions, AU- 
wide smooth variations, rather than sharp frontiers in the disk. 
Close to the sublimation temperature of the main dust elements, 
the disk mid-plane temperature becomes constant. This could re¬ 
sult in a less sudden and strong change of phase than initially 
thought, involving partial sublimation. Given an element in Ta¬ 
ble Q] and if we focus on the disk mid-plane region for which 
the temperature corresponds to that element sublimation temper¬ 
ature within a few Kelvin, the gas-to-dust ratio of that element 
would vary from ~ 0 at the outer edge of the temperature plateau 
to ~ 1 at the inner edge (the rest of the material is unaffected at 
these temperatures). A planetesimal crossing such a region in¬ 
ward might then be able to melt layer by layer, with possible 
consequences on its mineral petrology in the case of a future 
recondensation. Another way to model the partial sublimation 
would be to consider a vertical gradient of gas-to-dust ratio for 
the element in question: the element would be sublimated from 
the mid-plane up to a given altitude that would vary with the 
radial distance in the plateau. The farther into the plateau, the 
higher that limit, and the greater the proportion of the material 
of the column that is sublimated. Such a refinement would affect 
the dust vertical transport as studied by ICies 3 (1201(1 for ex¬ 
ample. Possible combinations of vertical and radial gradients of 
gas-to-dust ratio could also reflect the physics of the temperature 
plateaux. 

From these density and temperature radial profiles, we esti¬ 
mated the resonant torques that potential planetary cores would 
exert on the disk. From the total torque sign (resulting from the 
Lindblad and corotation resonances), we derived the locations of 
planetary traps and deserts and followed their own migration as 
the disk ages: after 1 Myr, four planetary traps were identified 
around 0.33, 0.60, 1.7, and 8.0 AU, with three planetary deserts 
around 0.41, 0.8 5 and 2.2 AU 

As stated bv iMasset et alJ ( 2006 ). density and temperature 
discontinuities appear to be the key characteristics leading to 
the generation of planetary traps: temperature troughs, density 
bumps, changes of phases, and heat transitions clearly affect the 
potential planet migration. Previous works indicated that den¬ 
sity and temperature irregularities can slow down type I migra¬ 
tion to allow planets to grow in a reasonable time. Our evolution 
code shows that such structures naturally arise during the disk 
evolution when the coupling between the disk dynamics, ther¬ 
modynamics, and geometry is taken into account: planets can 
be trapped at the main dust component sublimation lines. Ad¬ 
ditional simulations are still required to estimate by how much 
planetary embryos can be slowed down by being trapped, and 


whether they can actually migrate along with the planetary traps. 
It appears necessary to synchronize the planet and trap migra¬ 
tion rates, and this study requires a more thorough calculation of 
the planet feedback on the disk. Other sources of discontinuities, 
such as deadzones, may affect these migration torques and create 
more planetary traps or deserts. This will be treated in a separate 
paper. 
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